#!/usr/bin/perl -w
##This script run wiggler to convert bam file into visualized wig file, and zip them into binary wiggler (output_prefix.bw) file. 
##need to have fragment length estimation file from HOMER


## author: Yaping Liu  lyping1986@gmail.com 
## time: 2013-5-10

#Usage:  perl wrap_wiggler.pl [option] input.bam homer_dir output_prefix

use strict;
use Getopt::Long;

sub usage {

    print STDERR "\nUsage:\n";
    print STDERR "perl wrap_wiggler.pl [option] input1.bam input2.bam ... \n\n";
    
    print STDERR "  [Options]:\n\n";
	print STDERR "  --frag_len NUM : fragment length to extend the reads, allow multiple number, but the order should be the same as input bam or bed files. (Default: --frag_len 1, no extension)\n\n";
	print STDERR "  --prefix STRING : output file prefix(Default: --prefix sample)\n";
    
    exit(1);
}

my $WIGGLER = "/export/uec-gs1/laird/users/yaping/software/external/align2rawsignal/bin/align2rawsignal";
my $MAPPABILITY = "/export/uec-gs1/laird/users/yaping/data/genome_data/mappability/hg19/globalmap_k20tok54";
my $REFDIR = "/export/uec-gs1/laird/users/yaping/data/genome_data/genome/hg19/ucsc/maleByChr";
my $WIG2BW = "/home/uec-00/shared/production/software/UCSC_Browser_Tools/default/bedGraphToBigWig";
my $CHROMSIZE="/home/uec-00/shared/production/software/UCSC_Browser_Tools/default/hg19.chrom.sizes";
my $SORTBED = "perl /home/rcf-40/yapingli/storage/code/mytools/perl/sortByRefAndCor.pl";
my $REFINDEX="/export/uec-gs1/laird/users/yaping/data/genome_data/genome/hg19/hg19_rCRSchrm.fa.fai";

my @fragment_len = ();
my $prefix = "sample";


GetOptions( "frag_len=i" => \@fragment_len,
			"prefix=s" => \$prefix,
);

my @inputs = @ARGV;

usage() if ( scalar(@ARGV) == 0 );

if($#fragment_len == 0 && $fragment_len[0] eq ""){

	$fragment_len[0]=1;
}

my $wig_output=$prefix.".bedGraph";
##wiggler could only works correctly in BASH environment
my $cmd = "$WIGGLER -of=bg -n=5 -mm=7 ";
for (my $i=0; $i <= $#inputs; $i++){
	my $input = $inputs[$i];
	if($inputs[$i] =~ /\.bed$/){
		$input =~ s/\.bed$/.tagalign/;
		my $tmp = "ln -s $inputs[$i] $input";
		system($tmp);
	}
	
	my $frag = $fragment_len[$i];
	$cmd .= "-i=$input -l=$frag ";
}
$cmd .= "-s=$REFDIR -u=$MAPPABILITY -o=$wig_output\n";

print STDERR "$cmd\n";
system($cmd);

##sort the output bedGraph file
my $sorted_output=$prefix.".sort.bedGraph";
$cmd="$SORTBED --k 1 --c 2 --tmp ./ $wig_output $REFINDEX > $sorted_output\n";
print STDERR "$cmd\n";
if(system($cmd)==0 && -e $sorted_output){
	unlink($wig_output);
}

##convert bedGraph to big wig format
my $bw_output=$prefix.".bw";
$cmd="$WIG2BW $sorted_output $CHROMSIZE $bw_output\n";
print STDERR "$cmd\n";
if(system($cmd)==0 && -e $bw_output){
	unlink($sorted_output);
}


##wig to bigwig cost 1.5 times memory than the uncompressed .wig file.
#my $bw_output=$prefix.=".bw";
#$cmd="$WIG2BW $wig_output $CHROMSIZE $bw_output\n";
#print STDERR "$cmd\n";
#if(-e $bw_output && system($cmd)==0 ){
	#unlink($wig_output);
#}

